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Abstract. We describe fitting methods developed to analyze fluctuations in the Lyman-a 
forest and measure the parameters of baryon acoustic oscillations (BAO). We apply our 
methods to BOSS Data Release 9. Our method is based on models of the three-dimensional 
correlation function in physical coordinate space, and includes the effects of redshift-space 
distortions, anisotropic non-linear broadening, and broadband distortions. We allow for inde- 
pendent scale factors along and perpendicular to the line of sight to minimize the dependence 
on our assumed fiducial cosmology and to obtain separate measurements of the BAO angular 
and relative velocity scales. Our fitting software and the input files needed to reproduce our 
main BOSS Data Release 9 results are publicly available. 
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1 Introduction 

The surprising discovery [1, 2] of accelerating expansion in the current universe reveals that 
either some form of dark energy is driving the expansion or else that our theory of gravity 
is incomplete on the largest scales (see reference [3] for a recent review). The length scale of 
baryon acoustic oscillations (BAO) imprinted at the moment of baryon-photon decoupling 
provides a standard ruler that has been well measured in the temperature anisotropies of 
the cosmic microwave background [4, 5] (CMB) and in the number-density fluctuations of 
galaxies at z < 1 [6-24]. The Baryon Oscillation Spectroscopic Survey [25] (BOSS) of the 
third generation of the Sloan Digital Sky Survey [26] (SDSS-III) recently announced its 
ninth data release [27] (DR9), including an unprecedented number of high-redshift quasar 
spectra. The BAO feature is imprinted in the correlated fluctuations of intergalactic Lyman- 
a absorption of the light from these quasars, enabling us to measure the BAO standard ruler 
at redshifts z ~ 2.4, during the predicted era of matter-dominated deceleration. 

An optimal extraction of the BAO signal from a large sample of Lyman-a forest pix- 
els requires new fitting techniques beyond those previously developed to fit the large-scale 
clustering of galaxies [28]. For example, the non-uniform sampling of the absorption field, 
large expected redshift-space distortions, and relatively large depth suggest that estimates of 
the correlation function or power spectrum should be fit to three-dimensional models (rather 
than low-order multipoles at fixed redshift) expressed directly in terms of the physical ab- 
sorption wavelengths and angular separations between lines of sight. We describe here new 
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fitting techniques developed specifically for a near-optimal analysis of the DR9 correlation- 
function estimates described in a companion paper [29], and highlight the challenges and 
lessons learned. We are also making our fit input files and fitting code publicly available as 
a companion to this paper, so that readers may reproduce the main results presented here 
and in ref. [29]. 

The outline of our paper is as follows. In Section 2 we define our constituent models 
for linear theory with redshift-space distortions, non-linear effects, and redshift evolution. 
We also describe our parametrizations of possible deviations between our assumed fiducial 
cosmology 1 and the cosmology preferred by our data, and of broadband distortion of our 
correlation-function estimate introduced by analysis systematics. Finally, we introduce two 
data-reduction techniques using interpolated models. In Section 3, we describe the DR9 fit- 
ting inputs which consist of N = 1512 correlation- function estimates on a three-dimensional 
physical coordinate grid, accompanied by an initial estimate of their covariance in a block 
diagonal form that reduces the number of non-zero elements from N(N + l)/2 ~ 1, 114K to 
~ 64K. Next, we describe a novel method for internally validating and refining our initial 
covariance estimate, necessitated by a lack of available simulated mock statistics. In Sec- 
tion 4, we present our results on the fitting method and what it teaches us about the BOSS 
DR9 Lyman-a forest dataset. We discuss the expected parameter sensitivities and relative 
contributions from different regions of the three-dimensional separation space, and present 
model-independent data reduction results. Cosmological fitting results for DR9 are presented 
in the companion paper ref. [29]. We conclude in Section 5 with a discussion of the main 
lessons learned and our plans for future development. Appendix A provides details on public 
access to our fit inputs and fitting code. 

2 Models and Parameters 

We model a measurement of the correlation function 

t(r lfi ,z) = <*(si)<5(s 2 )) - (S( Sl ))(S(s 2 )) (2.1) 

where the ensemble averages are taken over realizations of a (possibly biased) tracer 5(s) of 
the large-scale distribution of matter in redshift space s, with (r, fj,, z) defined via 2 

1 f z dz' 
|s 2 -si|=r , |s 2 | - |si| = n-r , - |si + s 2 | = c / - (2.2) 

2 Jo H ?id\ z ) 

for some fiducial cosmology with Hubble function H^(z). The model combines a cosmological 
prediction ( cosmo with a parametrization of possible multiplicative and additive broadband 
distortions introduced by the analysis method. 

2.1 Physical Coordinates 

The physical coordinates for a pair of pixels measured in the absorption spectra of two 
quasars are the separation angle A8ij between the quasar lines of sight (A9{j = if the 
pixels are taken from the same quasar's spectrum) and the observed absorption wavelengths 

We assume a flat ACDM universe with ft A = 0.73, h = 0.7, fl h h 2 = 0.0227, and n s = 0.97 throughout. 
2 We use the notation fi = z • f and fik = z ■ k. 
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Aj and Xj. We convert the observed wavelengths to a relative velocity for the absorption 
systems 3 [30] 

Avij = c log {Xj/Xi) (2.3) 

and an average absorption redshift 

Zij = ^5 _ l (2.4) 

where X a ~ 1216 A is the rest wavelength of the Lyman-a transition which determines pixel 
redshifts Zi = Xi/X a — 1. We calculate the corresponding co- moving separations along (rn) 
and perpendicular to (rj_) the line of sight as 



1 + z 



r\\ = D CM { Zj ) - D CM { Zi ) = - - l \ ■ A Vij 1 + O (Av^/cY (2.5) 

tlM{Zij) L J 

f Zi i dz' 

r± = D AM ( Zij ) ■ Adij = c / — — —— ■ AOij (2.6) 

where Dq^(z) and -D^fid^) are the co-moving line of sight and angular distance functions, 
respectively, for the assumed fidicual cosmology H&&{z). We introduce separate scale factors 
aii and a±_ in Section 2.2.4 to allow for small discrepancies between the true cosmology 
and our assumed fiducial model. Fig 1 shows that the BOSS blue camera wavelength limit 
of A > 3600 A limits pixel pairs contributing to a BAO peak feature near 110 Mpc/h to 
AOij ^$ 100 arcmin. Similarly, the BOSS redshift coverage Zij < 3.25 limits the BAO peak to 
Avij/c = log(A2/Ai) < 0.04 and observed wavelength differences |Aj — Aj| < 200 A. 

2.2 Cosmological Models 

We build the cosmological model starting from an isotropic linear power spectrum prediction 
P(k,zo) at some reference redshift z$, then embed this prediction in redshift space (we use 
tildes to denote linear-theory predictions without any redshift space distortions). In the 
general case of a plane-parallel redshift-space distortion r — > (r, (j) we have [31]: 

Ccosmo^, H,Zq) = Li{p)£l,casma{T,Zo) (2.7) 

£even 

with 



i poo 



&,cosmo(r>*o) = J k 2 j e (kr)P e (k,z )dk (2.8) 

where Li is the Legendre polynomial, ji is the spherical Bessel function, and Pi(k,zo) are 
the multipoles of the redshift-distorted power spectrum P(k, Hk, zq) with = z • k. 

2£ + 1 f +1 

P e (k,z ) = — - — J P(k,iJ,k,za) Li(nk)dfJ,k • (2-9) 

Specializing to linear theory and the distant observer approximation [32], the infinite series 
of eqn. (2.7) is truncated at £ = 4, with 

P e (k, z ) = b 2 (z )C e ((3(z ))P(k, zq) (2.10) 



3 This definition is not identical to the Doppler velocity that an observer at one absorber would measure 
for the other absorber, but does agree to second order in the wavelength ratio. 
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Figure 1. Physical coordinates for pixel pairs with observed Lyman-a absorption wavelengths A.; < 
Xj. Grid lines of Zij (vertical blue, values left to right are 2.25, 2.75, 3.25), Avij/c (horizontal red, 
values bottom to top cover 0.001-0.049 with 0.002 spacing, with additional contours at 0, 0.059, and 
0.083) represent the nominal sampling grid used in a fit. The shaded gray region shows the pixel 
pairs contributing to a typical BAO fit, bounded by Ai > 3600 A , r» < 170 Mpc/h, and Zij < 3.25. 
Contours of Ad (thick black curves) at which the 3D separation is 110 Mpc/h (values from bottom 
left corner out are 100, 80, 60, 40, 20, arcmins) identify pixel pairs contributing to the BAO peak 
region at different angular separations. 



and 



Q03) 



2£ + l 



(l + /3^l) 2 L e (/j,k) d/j, k = < 



l + |/3 + l/3 2 
|/3 + |/3 2 

8 o2 



-2-/3 



(2-11) 



where b(z) and j3{z) are the tracer bias and redshift-space distortion parameter at redshift 
z, respectively. We can therefore write 



£cosmo(r,/x,zo) = b 2 (z ) Ce(P(zb))L£((j,)£e >coamo (r, 



(2.121 



£=0,2,4 

in terms of the undistorted linear-theory multipoles 

■£ />oo 



&,c»smo(r,Zo) = ^ J k 2 j £ (kr) P(k,z )dk 



(2.13) 



Figure 2 shows examples of the linear models we use for Lyman-a fitting in this paper, 
calculated with zq = 2.25. Note that the undistorted multipoles ^ are not independent since 
they derive from the same underlying power spectrum via eqn. (2.13). Explicitly, we find 
that £2(0 arid £4(7*) can be calculated directly from £o(r') specified on an interval ro < r' < r 
via (we have dropped zq here for clarity): 



= £o(r) + 



6(»"o) - £o(ro) 



£4(7-0) - £0(7-0) 



£o(r')r /2 dr' 



(2.14) 



r' A dr' 



-A- 




Figure 2. Cosmological linear models calculated for zq = 2.25 and assuming a flat universe with 
f^A = 0.73, h = 0.7, il^h 2 = 0.0227, and n s — 0.97. Panels show the (fc-weighted) power spectrum 
(top-left) and the (r 2 -weighted) correlation function monopole (top-right), quadrupole (bottom- left), 
and hexadecapole (bottom-right). Curves are calculated with CAMB [33] (thick, red) and using ref. [34] 
(light, blue) with solid curves showing the full cosmological model and dotted (dashed) curves showing 
the corresponding CAMB "sideband" ( "no- wiggles" of ref. [34] ) smooth model. 



2.2.1 Peak Decomposition 

We expect a Lyman-a analysis to distort the measured broadband shape of the correlation 
function, so our goal is to only use information from the localized peak near r — 110 Mpc/h 
to measure BAO parameters. To achieve this goal, it is useful to decompose ^ cosmo into 
separate "peak" and "smooth" (or "no- wiggles" ) components. Following ref. [34], we can 
isolate the oscillations in a multiplicative term of the effective baryon transfer function Tb(k) 
(see equation (16) in ref. [34]), which leads to a rather complicated decomposition for the 
correlation function. Instead, we adopt the unphysical but more tractable decomposition 

£e,cosmo(r, Zq) = ^smooth (r, Zq) + 6?,peak(j, Zq) (2.15) 

with the understanding that this model is only valid when using parameter values that are 
sufficiently close to their nominal values to recover a physically plausible ^ iCOSmo . 

Figure 2 shows the smooth P(k, zq) and corresponding correlation multipoles suggested 
in ref. [34]: the oscillations are effectively removed, but the corresponding "peak" feature, 
defined as the difference between the full and "no-wiggles" models, is not well localized in 
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any of the multipoles, with deviations from the dashed curves extending far from the peak. 
To remedy this problem, we construct an alternate peak model that is explicitly localized, 
but only for a single linear combination of multipoles (or, equivalently, a single value of /i) 
because the integrals in eqn. (2.14) effectively spread the peak to all scales above and/or below 
(depending on the choice of Tq) for other linear combinations. Our CAMB "sideband" smooth 
model is constructed as follows: we first isolate a localized peak in the CAMB prediction [33] 
of Figure 2 by simultaneously fitting the regions 50-86 and 150-190 Mpc/h of the monopole 
to the form 

+1 

£o,fit(r,z ) = E c i rJ ■ ( 2 - 16 ) 
i=-3 

Next, we replace £o,cosmo with £o,fit in the region 86-150 Mpc/h to obtain £o,smooth- We then 
calculate |2, P eak and | 4 , P eak using eqn. (2.14) with &,peak(r-o) = at r = 0. Finally, we 
calculate 

^.smooth (r, Z ) = & lC osmo(»", Z ) - & lPea k(r, Z ) ■ ( 2 -17) 

The resulting CAMB "sideband" smooth model is shown in Figure 2. 

By construction, our CAMB peak multipoles are exactly zero below 86 Mpc/h, but the 
peak spreads to large scales for I = 2,4, as required by eqn. (2.14). Note that although it 
would be possible to specify independent localized peaks for each multipole, this is unphysical 
(within the framework of ref. [32]) and does not, in general, reproduce the evolution of the 
peak shape in £(r, h,zq) with fj, implied by the constraint of a single underlying P(k,zo), 
as illustrated in Figure 3. For example, the constrained evolution leads to a percent-level 
shift in the peak position as a function of /i shown in Figure 4 for the two models used here. 
Note that this shift is essentially the same for both models, despite the rather different peak 
shapes shown in Figure 3. 

0.0015 

- 0.0010 

In* 

^ 0.0005 
0.0000 



50 



Figure 3. Cosmological peak models calculated using the CAMB "sideband" method described in 
the text (left, red) and the "no-wiggles" method of ref. [34] (right, blue) described in the text. Curves 
show fi = 0.4 (solid), 0.7 (dashed), and 1.0 (dotted). There is no r 2 weighting applied here. 

2.2.2 Nonlinear Effects 

The expected effects of non-linear structure growth on the BAO feature can be modeled with 
an anisotropic Gaussian roll-off of the linear power spectrum [35] : 

PNL(k, fi k , Zq) = exp(-k 2 E 2 (fx k )/2) ■ P(k, z ) (2.18) 



0.0010 




200 



r (Mpc/h) r (Mpc/h) 
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Figure 4. Evolution of the peak position with /i for the cosmological peak models shown in Figure 3. 
Curves show f3 = 1.0 (solid), 1.4 (dashed), and 1.8 (dotted). Fractional shifts are measured relative 
to the position of the monopole peak for each model. 



where 



S 2 ( / x fe ) = / x^ + (l-^)Ei. (2.19) 



In general, this approach breaks the decomposition, eqn. (2.10), of Pi(k,zo) into separate j3- 
and fc-dependent factors and requires that the integrals of eqns. (2.8)-(2.9) be re-evaluated 
for each value of f3. However, since we expect E ~ 5 Mpc/h, compared with an expected 
peak full- width half- maximum of ~ 25 Mpc/h, we can approximate for f3 ~ (3q: 



with 
and 

ft(P) = 



Pe,NL(k, z ) ~ exp(-A: 2 S £ 2 (/3o)/2) • P e (k, z ) (2.20) 
S 2 (/3) = f e (J3) ■ Ejj + (1 - //(/3)) • Si (2.21) 

f-l f4 i 1 + P^lf L dVk) dfi k 



35+42ff+15/3 2 
105+70/3+21/3 2 ' 
7+12/3+5/3 2 

14/3+6/3 2 
15 , 2 f = 4 





2 . (2.22) 



This approximation effectively models anisotropic broadening using different amounts of 
isotropic broadening for each multipole. The resulting correlation function multipoles £e,NLy 
calculated with eqn. (2.8), can then be substituted in eqn. (2.12). 

Taking a fiducial value of /3q = 1.4, we calculate /o = 0.505, /2 = 1.07, and = 
2.79. At a redshift z = 2.4, we expect Ey ~ 6.41 Mpc/h and E ± ~ 3.26 Mpc/h, so that 
So = 5.10 Mpc/h, E2 = 6.58 Mpc/h, and E4 = 9.79 Mpc/h. Figure 5 compares the 
resulting approximate models with exact calculations for a range of (5 values. Note that we 
are neglecting the non-zero even multipoles I = 6,8,... that are introduced by anisotropic 
non-linear broadening with this approximation, but Fig 5 shows that most of the BAO peak 
feature has already been smoothed out by I = 4. We also neglect any redshift evolution of 
the non-linear broadening parameters, since this evolution would be a second order effect 
on what is already a small correction to the linear theory for the purposes of measuring the 
BAO feature. 

In order to achieve a self-consistent decomposition of each correlation-function multipole 
into peak + smooth components with non-linear broadening effects included, we first broaden 
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Figure 5. Effects of anisotropic non-linear broadening implemented with eqn. (2.18) and applied 
to z = 2.25 linear CAMB predictions [33] using S|| = 6.41 Mpc/h and E ± = 3.26 Mpc/h. Curves 
show no broadening (thick red, same as curves in Figure 2), isotropic broadening (dotted red) by 
(S| + E^) 1//2 /2 = 5.09 Mpc/h, the approximate anisotropic model described in the text (dashed 
blue) with j3o = 1.4, and the envelope of full anisotropic calculations (light blue shaded) for /3 = 0.5- 
2.5. Left-hand panels show the fc-weighted multipoles Pe,NL(k, Zq) with b 2 Ce(f3) divided out, for I = 
(top), t = 2 (middle), and I = 4 (bottom). Right-hand panels show the corresponding r 2 -weighted 
correlation function multipoles £^(r, Zq). 

each ^£ tPea k(r, Zq) by E^ to obtain ^, pe ak,AfL(^, zq). Broadening by E in this context implies 

P(k) -> P(Jfe) exp(-A: 2 E 2 /2) (2.23) 

which transforms to 

£o(r)-> / da- [G(r-s,E)-G(r + s, £)]&(*) , (2-24) 
jo r 
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where G is the normalized one- dimensional Gaussian 

Note that when r > S, we recover the expected convolution 

/>oo 

£o(r) -> / da G(r- a, (2.26) 
J o 

to a good approximation. Figure 6 shows the broadened multipole peaks based on our 
linear CAMB templates. To derive the corresponding smooth templates, we either use the 
linear smooth templates, or else we calculate non-linear smooth templates by subtracting the 
broadened peak templates from the broadened multipoles: 

&,smooth,iVL = £t,NL ~ &,peak,iV£ • (2.27) 

In the first case, we are only applying non-linear broadening to the peak feature, which 
introduces an unphysical distinction between peak and smooth components. In the second 
case we are broadening the full correlation function, which removes the distinction but applies 
unphysical filtering of small scale structure. A better model is probably somewhere in between 
so we consider both alternatives, and expect that any extracted peak parameters will not 
depend on this choice. Figure 7 compares both approaches. 




r(Mpc/h) r ± (Mpc/h) 



Figure 6. Effects of anisotropic non-linear broadening on peak templates, implemented using the 
approximations described in the text and applied to zq — 2.25 linear CAMB predictions [33]. Solid 
(dashed) curves show peak templates with (without) non-linear effects. From top to bottom, curves in 
the left-hand panel show the I = (red), 2 (blue), and 4 (green) correlation multipoles, with b 2 Ce(f3) 
divided out, without any r- weighting and scaled by 10 4 . All templates are identically zero for r < 60 
Mpc/h, by construction. The right-hand panel shows equally spaced contours of £ P eak(?"j fi) calculated 
with j3 = 1.4, with the outer contour corresponding to zero correlation. 
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Figure 7. Comparison of two different schemes for incorporating anisotropic non-linear effects into 
the final templates used for fitting. Panels show the (fc-weighted) power spectrum (top-left) and the 
(r 2 -weighted) correlation function monopole (top- right), quadrupole (bottom-left), and hexadecapole 
(bottom-right). Red curves are calculated with non-linear effects applied to the BAO peak feature 
only. Blue curves are calculated with non-linear broadening applied at all scales. Both schemes are 
derived from zq = 2.25 linear CAMB predictions [33] and use the same broadened peaks shown in 
Figure 6. The combined peak + smooth (smooth only) templates are represented with solid (dashed) 
curves. 



2.2.3 Redshift Evolution 

In general, we model the redshift evolution of a parameter p(z) for z near zq in terms of two 
parameters po and j p via 



p( z ) = Po [ — ) • ( 2 - 28 ) 
\1 + z J 

We apply this evolution to the parameters b 2 (z) and (3(z), introduced above, and to the BAO 

scale parameters ai SO (z), ot\\(z) and a±(z) introduced below. Given a covariance matrix for 
the parameters po and 7 P , 

P<T ^\ (2.29) 
POqo 1 a* J 



Cp 



7 

the variance of p{z) with z ~ z$ is given by 

al {z) = J-C v -J l (2.30) 

where 

V dp d-f p J 



is the Jacobian. The error on p(z) is smallest at 

log (^-) =-b-pa + V& 2 - (1 - P 2 W (2-32) 
\ l + Zo J 



with 

(TO ,1 

a = , b = — 

2.2.4 Scale Factors 



(2.33) 



When fitting our model to data, we allow for an overall relative normalization factor a pea k — 1 
as well as a general coordinate transform, r' = r'(r,fi,z) and p! = fj,{r, fj,, z), that allows for 
possible small differences 4 between the fiducial and actual cosmologies: 

^cosmo('*! Pt z) y Apeak " [^{T > M ) Z) ^smooth (j" i A* > ^)] "I - Csmooth(^ 

f ,f/ ,z) . (2.34) 

We consider two options for the cosmological broadband £smooth : we either apply the same 
transform as for the peak (r" = r'), or else we keep it fixed (r" = r) so that only the 
peak is transformed. Neither approach is exact when considering variations of the cosmo- 
logical parameters around our fidicual model, but we find that decoupling the peak from 
the cosmological broadband (r" = r) better localizes the separations contributing to a BAO 
measurement to the peak region (compare Figs. 16 and 19), so is preferred when broadband 
distortion is not fully under control. The corresponding /c-space transforms are defined by 

k'(k,n k ,z) -r'(r,fj,,z) = k-r , f/ k (k, fi k , z) ■ fj,(r, fx, z) = fi k ■ n . (2.35) 

Including a /^-dependence in the coordinate transform also enables us to study the con- 
straining power of our data separately along and transverse to the line of sight. Similarly, 
a ^-dependence allows us to determine the redshift at which our results are best measured 
using eqn. (2.32). 

For our baseline isotropic model, we use 

r[ m {r,n,z) = a iso {z) ■ r 

/ / \ ( 2 - 36 ) 
Mi so ( r >M = A* • 

For our baseline anisotropic model, we decouple the line-of-sight (rii — y cnurn) and transverse 
(rj_ — y a±r±) scales using 

r'^ir, n, z) = a ani (p,z)-r 

, . , <*„(*) (2.37) 

OanilAi, Z) 

with 

aaniGu, z) = ^a\{z)^ + a\{z){l-^) . (2.38) 

When the parameters of this anisotropic model are determined by measuring the physical 
coordinate separation scales Avbao(-z) and A#bao(z) corresponding to the comoving BAO 
scale re ao a t some redshift z, related by 

Av BA o(z) ~ r BAO H(z)/(1 + z) (2.39) 
A^bao(^) = r BAO /D A (z) , (2.40) 



4 In case there is evidence for large differences, the analysis should be repeated with a fidicual cosmology 
that better matches the data. 
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then the expected best fit values of an and a± satisfy 

1 + 2 

"II • u i n At;baoO) - ^BAO.ad (2.41) 

a± ■ D a ,m( z ) A °bao(z) = r B AO,fid (2.42) 

where r B Ao,fid(^) is the comoving BAO scale predicted by the assumed fiducial cosmology. 
Combining these results, the best-fit anisotropic scale factors measure 

, s r-BAO,fid H M (z) r BA o,fid D A (z) 

«||(z) = TjTX~ ' a -L z ) = n • 2 ' 43 

11 r B AO H[z) r B AO D AM (z) 

The distributions of these ratios predicted at z = 2.4 by WMAP9 [4] observations of the 
CMB are shown in Figure 8. 




Figure 8. Predictions at z = 2.4 for the anisotropic scale factors (left) and a an i(/j,) (right) from 
WMAP9 data [5]. Filled blue contours show ACDM predictions, while solid (dashed) curves show 
the effects of allowing the curvature ttk (dark energy equation of state parameter wq) to vary in the 
model. The gray filled contours show the combined effects of varying both ilk and wq. Contours 
in the left-hand plot enclose 68% and 95% while the bands in the right-hand plot are ±1 standard 
deviation. The fiducial cosmology ot|| = a± = a(fi) = 1 is indicated by dotted lines. 



We can evaluate a transformed correlation function £' in the original coordinates (r, /x, z) 
using eqn. (2.12) (with z-dependencies omitted for clarity) 

cosmo 

(r'(r, ii),n'{r, fi)) 
= b 2 £ C^)L^\r^))^ cosmo {r'(r^)). ( 2 - 44 ) 

£=0,2,4 

The corresponding multipoles in the original coordinates are then given by 
of _i_ i r+ l 

£(r) = — ^- J t>{r,n)L^)dv (2.45) 



2 

2^ + 1 



b 2 / h,co Sm o{r'{r^))L^(r^))-L^)d^. (2.46) 

£'=0,2,4 
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In the case of isotropic distortion, we can simply replace r with a; so (z) • r in the multipoles 
obtained with eqn. (2.13). However, anisotropic distortions mix multipoles and lead to more 
complicated expressions. For small e = (an — aj_)/(ai| +ot±) and arbitrary a = (an +aj_)/2, 
we have 

an=a(l + e) , a± = a(l - e) (2.47) 



and corresponding transforms 

4(r,/i) = ar [1 + 2e(2/i 2 - 1) + e 2 ] 1/2 ~ ar [l + (2,u 2 - 1) • 



£'=0,2,4 



where 



(2.48) 



/ u / ani (r, / x) = / i(l + e)[l + 2e(2^ 2 -l) + e 2 ] 1 ~ M [l + 2(1 - /x 2 ) • e] . 
Using these approximations, and 

6,cosmo(? ,/ (^, m)) — C£,cosmo(ar) + e • ar • d r &/ cosmo (ar) • (2fi 2 - 1) , (2.49) 
we obtain distorted multipoles: 

£e( r ) - 1)2 ^2 C £'(^ {^,cosmo(H • A e ^(e) + e • ar ■ d r £ e > iCOSmo (ar) ■ B^/j (2.50) 



2/ + 1 

^'(e) = y i Ma* + 2e/x(l - ti 2 ))L e (fj,)dfi (2.51) 

Btjt, = J ^ (2fi 2 - l)L e (fi)Lt(n)dn . (2.52) 

We find the following non-zero coefficients, to first order in e: 

14 4 
A),o = 1 , B a = — - , A _ 2 = -e , B 02 = — (2.53) 

3 5 15 

4 4 1 40 8 

#2,0 = 3 , ^2,2 = 1 + 7 e > ^2,2 = ^ , A 2A = —e , £ 2 , 4 = — (2.54) 

48 24 40 1 

A4,2 = -gg€ , ^4,2 = ^ , ^4,4 = l + ^e , B 4,4=^ ( 2 '55) 

80 20 , 

6 ' 4 = _ 33 e ' M = 33' (2 ' 56) 

Note that an £ = 6 term must be included in order to match eqn. (2.46) to first order in 
e, but is numerically negligible. On the other hand, second-derivative terms in the Taylor 
expansion 2.49 are formally 0(e 2 ) but are numerically significant when 

e > 2| ^1 (OT)I , (2.57) 
ar\d*£e(ar)\ 

which is guaranteed to occur at the peak of the BAO feature and at the transition to a 
rising smooth broadband below the peak. Figure 9 shows examples of anisotropic coordinate 
transforms based on eqns. (2.44) and (2.50). Note that eqn. (2.50) is only exact when fi 2 = 1/2 
(rii = rj_) and grows less accurate when moving away from the diagonal, while eqn. (2.44) is 
exact and therefore preferred when fitting for the anisotropic scale parameters. 
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r(Mpc/h) r ± (Mpc/h) 

Figure 9. Examples of anisotropic coordinate transforms applied to the correlation function, using 
j3 = 1.4. The top row is calculated with a — 0.95, e = 0.05 and the bottom row with a = 1.05, e 

0.0"). both chosen to match the fiducial cosmology at fi = 1, where the expected sensitivity is largest 
for j3 — 1.4. Left-hand plots show the normalized distorted correlation multipoles 10 4 ^(r)/(6 2 C^ (/?)) 
calculated with eqn. (2.50) (solid curves) or without any coordinate scaling (dashed curves). From 
top to bottom, curves show I = (red), 2 (blue), and 4 (green). Right-hand plots show equally spaced 
contours of £(r,fA,Zo) with the outer contours corresponding to zero correlation. Thick solid curves 
are calculated with eqn. (2.44) and dashed curves are calculated without any coordinate scaling. 
Thin magenta curves show the results of combining the transformed multipoles (solid curves) of the 
left-hand plots using eqn. (2.8), and demonstrate the level of accuracy provided by the first-order 
approximations described in the text. 
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2.3 Broadband Distortion Models 

A Lyman-a analysis measures the correlation function 



to - <« - {km = {F %^ Srr 1 (2 - 58) 



for the overdensity proxy 

. Ft - F(Xi) 

Oi = = 

in pixels i measured at rest absorption wavelengths Aj, where 

Si /sky,i 

* = ~f ~~ 

J cont,j 



(2.59) 



(2.60) 



is the transmitted flux fraction, relative to the sky flux level / s k y ,i and normalized to an 
assumed quasar continuum flux level / C ont,i> an d F(Xi) is an assumed mean transmission 
fraction at Aj. 

To investigate possible sources of broadband distortion, write 

fi = fi + £i 

/sky.i = /sky.i + £sky,i + s(Aj) (2-61) 
/cont,i = /cont,i ~i~ Cj(d) 

where tilde quantities are true values, e, and e s k y ,i are (zero-mean) noise sources, s(Aj) ac- 
counts for any wavelength-dependent residual sky-subtraction bias (as observed in DR9), and 
Cj(d) describes the continuum modeling error that, in general, depends on the full vector of 
pixel measurements 

d = {*j,fj}?=r ela • (2-62) 
Combining eqns. (2.60) and (2.61), we find: 

p fi /sky,i 'S(Aj) ~\~ £j ^sky,i 

/conM+Ci(d) (2.63) 

~Fi[l-Si-Ci + Ek] 

with 

p.= fi~ iW,* i Si = j( X i) ; Q = -^- , ^ = S ~ , (2.64) 

/cont,j /i /sky,i /cont,i /i fsky,i 

where we have assumed that Cj <C 1. Averaging over noise and cosmic realizations, we find: 



to — Tij 

where 



to (1 + Aij) + Bij (2.65) 



- _ (F l F 3 )-(F l )(F J ) 
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is the undistorted cosmological correlation function we seek to measure, rjj ~ 1 describes the 
effects of any discrepancy between the assumed and true mean transmitted flux fraction via 

- - (2.67) 



Aij is a multiplicative distortion of the true correlation function 

Aij = (5,5,) + (dCj) + (EiEj) + (SJiCj) + (5 J )(C i ) - (S l ) - (Sj) - (d) - (Cj) , (2.68) 

and Bij is an additive distortion 

B i:i = (5,5,) - (Si)(Sj) + (QCj) - (Ci)(Cj) + (EiEj) , (2.69) 

where we have used (Ei) = and assumed that 5, C, and E are mutually uncorrelated. 

Equation (2.65) demonstrates that, in general, the measured correlation function can 
be systematically distorted by multiplicative and additive effects. The actual correlation 
function used for fitting is a weighted sum of pixels in bins of (r,fj,,z). Any multiplicative 
distortion is due to a mismatch between the assumed F(X) and the true mean transmitted 
flux for pixels at the same wavelength A, whereas additive distortion is due to correlated 
continuum fit errors. 

We model broadband distortion using the following parametrization that combines the 
multiplicative (B m ) and additive (B a ) effects described above: 

/ 1 + z \ 7fc2 

£(r, /x, z) = £cosmo(r, /x, z)-[l + B m (r, /x, z)] + B a (r, /x, z) • — (2.70) 

\l + ZqJ 

with (B x represents either B m or B a ) 

/.(*) . ( r ^ . r.r.A.f 1 + z 

i= ^min J = Jmin n=ri m'm 

where, nominally, tq = 100 Mpc/h and zq = 2.25, and 



2 max Jmax ^ma 



Bx (r, ,, z) =y: e e c • (i - *«J • • I f^J ' (2 - 7l) 



i < 

1 t > 



(2.72) 



The integers i mm < i max , < j m \ n < j m£LX , and n mm < n max determine the number of free 
parameters, and we normally restrict j to even values. The t = 0,2,4 multipoles of the 
distorted correlation function are given by (arguments of (r, z) are suppressed for clarity) : 

£o = £o .cosmo (1 + -Bo.m) + £2 , cosmo , cosmo --B4,m + Bo 

5 9 

( 2 2 \ 

£2 = Co,cosmo - (1 + -02,m) + ^2,cosmo ' ( ^0,m + ^^2,m + -jBi,m j + 

£4,cosmo • f -jB2 t m + ^^-&4,m + 143 -^ 6 > m ^ ^ 2 > a (2.74) 
?4 = £o,cosmo • (1 + i?4,m) + ^2,cosmo • ^ — £?2,m + ^jB^ rn + J^-^6,mJ + 

C4,cosmo • f-Bo,m + ^^2,m + iQQi ^ 4 . m + 143 ^ 6,m ~*~ 2431 ^ 8 ' m ) ~*~ ^ A ' a ( 2 -75) 
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with distortion multipoles given by 

*max ^max / \ j / 1 i \ n 

E E -(rrl- 1 ) ■ < 2 - 76 > 

* — ^min ** ^min 

For the purposes of fitting DR9, we chose six different parameter configurations for broadband 
distortion, as summarized in Table 1 and plotted in Figure 10. 



Name 


i 


j 


n 


i 


B a 

j 


n 


BB1 








0,1,2 


0,2,4 





BB2 








-2,-1,0 


0,2,4 





BB3 








0,1,2 


0,2 


0,1 


BB4 


0,1,2 


0,2,4 











BB5 


0,1,2 


0,2 


0,1 








BB6 


0,1 


0,2,4 





0,1 


0,2,4 






Table 1. Broadband distortion models used to fit DR9. Models are labeled BB1-6 and columns show 
the range of indices in eqn. (2.71) used in the multiplicative (B m ) and additive (B a ) components. 
Dashes indicate that a component is not used. 




Figure 10. Best fits of the six distortion models BB1-6 described in Table 1 to DR9. Plots show 
the r 2 -weighted distortion £(r, fx, z) — Ccosmo^, (J., z) for z — 2.5. Models BB1-3 are on the top row 
of plots, left to right, and BB4-6 are on the bottom row. The faint negative impression of the BAO 
feature visible for BB4 and BBS reflects the fact that they model purely multiplicative distortions. 
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2.4 Interpolated Models 

For the purposes of visualization and studying non-cosmological signatures in data, it is use- 
ful to fit using simple interpolated models with minimal assumptions. We describe here two 
such models: the first treats the correlation multipoles as (possibly independent) arbitrary 
interpolations through fixed arbitrary values of the separation, and the second treats cor- 
rections to a smoothed power spectrum as (possibly independent) arbitrary interpolations 
through equally spaced wave numbers covering a limited band. 

For our first model, we pick n fixed values of the separation fa} and model the r 2 - 
weighted correlation multipoles as: 

r 2 &(r, z) = b\{z) I m , n (r, z ; fa}) (2.77) 

where I m ,n is an m-th order interpolation through the n points fa} evaluated at separation 
r and referenced to the redshift Zq. Comparing with eqn. (2.10), we write 

b 2 (z) = b 2 (z)C £ (f3(z)), (2.78) 

so that 

I m ,n(r, Zq) = r 2 £z, cosmo (r, Zq) (2.79) 

when £t(r,[J,,z) is purely cosmological. 

For our second model, we interpolate over a fixed band of fc-space and therefore require 
some assumption about the power outside this band in order to predict r-space correlations. 
We adopt the following form for the multipoles of the power spectrum: 

P e (k, z) = bj(z) [P sm ooth(k, Zq) + AP e (k, z )] (2.80) 

where -P S mooth(fc ) zq) is a smoothed fiducial power spectrum, without any BAO features, mod- 
ulated by functions APi(k, Zq) that are zero outside some interval (fciojMu), both referenced 
to zq. The correlation multipoles corresponding to eqn. (2.80) are obtained via eqn. (2.8): 

&(r, z) = b 2 (z) ie jSmoo th(r, zq) + A^(r, zq) (2.81) 

with 

Ae<(r,2b) = 2-2 / k 2 j e (kr)AP E {k,z )dk . (2.82) 

We parameterize the /c-weighted modulation as a linear combination of m-th order B-spline 
functions B m : 

n—m—2 / , , \ 

k-AP e (k,z )= hjB m ( " j ) (2.83) 

3=0 

with n uniformly spaced knots kj at: 

kj = h + j ■ Ak , Ak = khi - ho . (2.84) 

n — 1 

Using the property that B m (t) is only non-zero for < t < 1, we can write 

ACe(r,z Q ) = ^ Yl k,j E e,m(r;kj,Ak) (2.85) 
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where 



E lim (r- kj, Ak) = {m+ l)Ak [ k{t) j e (fc(i)r) B m (t)dt (2.86) 

Jo 

with 

k(t) = kj + (m + l)Ak ■ t (2.87) 

can be solved analytically in terms of trigonometric functions and the sine integral. Figure 11 
shows examples of Ei^ m for cubic (m = 3) B splines. 




Figure 11. Cubic (m = 3) basis spline contributions to the (/c-weighted) power multipoles kAPg(k) 
(top left) with their corresponding (r 2 -weighted) contributions -E^ m (r; kj,Ak) with ^ = 0,2,4 (clock- 
wise from top right) to the correlation function multipoles A^(r) with fc lo ,fc hi = (0.02,0.2) h/Mpc, 
n = 10, and j = (solid blue curves) and j = 5 (dashed red curves). 



3 Fitting Method 

We fit the parameterized models of £(r, z) described above to estimates = £,( r ijk, f^ijk, z k) 
of the correlation function specified on a 3D grid Avi ® A6j ® z k of physical coordinates (see 
Section 2.1), with 

r ijk = yJr^{Avi,z k ) + r\{Ae j , z k ) (3.1) 
Hijk = r\\{Avi,z k )/r ijk . (3.2) 

Our nominal coordinate grid consists of 28 unequally spaced points in At; covering < Av < 
0.083c (see Figs. 1 and 12), 18 equally spaced points in A6 spanning 5-175 arcminutes, and 
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three redshift values (2, 2.5, 3), for a total of N = 1512 grid points. Each set of correlation- 
function estimates is accompanied by estimated covariances 



for each physical coordinate pair (ijk) and (i'j'k'). In the following, we use the notation d to 
refer to a vector of N correlation function estimates £ijk, and write C for the corresponding 
N x N covariance matrix. See ref. [29] for details on the estimates d and C obtained from 
BOSS DR9 and accompanying simulated mock data that we use below. Appendix A provides 
instructions for downloading these estimates as well as the software necessary to reproduce 
the main results provided here and in ref. [29]. We define a standard chi-square in terms of 
a parameter vector and theory prediction fJ,(6) as 



3.1 SDSS-III Data Release 9 Inputs 

Figure 12 compares the data vectors for theory, simulated BOSS DR9 data, and actual BOSS 
DR9 data. The signal is visually obscured in the data by large-amplitude modes at fixed 
separation angles AO ~ r±. However, these modes have negligible impact on our ability 
to measure BAO parameters since they have very small weights (by construction) in the 
covariance matrix (see Figure 13). The value of log 10 \C\/N provides a rough measure of the 
overall (i.e., not BAO specific) signal to noise ratio of a sample. We find values of -8.5 and 
-8.3 for data and mocks, respectively, with N = 1512, which indicates that the errors in our 
simulated data are about 30% larger than in real data. The signal to noise ratio is largest at 
z = 2.5, with errors about 40% larger at z = 2 and three times larger at z = 3. 

Figure 13 shows the structure of our covariance estimates, which is essentially the same 
for simulation and data after accounting for the ~ 30% normalization difference. The first 36 
eigenmodes have artificially large eigenvalues (see Figure 13a) due to the template marginal- 
ization procedure described in ref. [29], which accounts for the largest-amplitude modes visible 
in Figure 12. We expect correlation function estimates to be uncorrelated between different 
separation angles and therefore impose this constraint on our estimated inverse covariance 
matrix reducing both C and C _1 to a block diagonal form consisting of 18 sub-matrices, Cu\ 
and C7j\ , with dimensions 84 x 84. Each of these submatrices has the same structure, shown 
in Figure 13c for C^, but a different normalization, \C^\ or \CZ\\ = | C ( j ) | 1 , depending on 
its separation AOj (solid curve in Figure 13b). Each of the three 28 x 28 sub- matrices on the 
diagonal, corresponding to a fixed separation AOj and redshift z^, also have the same struc- 
ture with normalizations shown as the dashed curves in Figure 13b. Values of the inverse 
covariance along the diagonal are roughly independent of Av% and only depend on A6j via 
the factor |Cq-)|, as plotted in Figure 13d where the decreasing values at large Avi correspond 
to increasing errors due to limited statistics. 

Our correlation function estimates are actually provided as a set of M quasi-independent 
data vectors d m and covariance matrices C m , corresponding to angular partition of the BOSS 
survey footprint into individual observing plates [29]. For fits to DR9, we have M = 817 
plates. We perform a weighted combination of observations, assumed to be independent, to 
obtain the final d and C used for fitting 



{iijk^i'j'k 1 ) - {iijk){ii'j'k') 



(3.3) 



(3.4) 



M 



M 





(3.5) 
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Figure 12. Slices in (r, fx) of the estimated r 2 -weighted correlation function for the fiducial cosmology 
(top row), simulated mocks (middle row), and BOSS data (bottom row). Grid lines show the actual 
sampling of used. From left to right, columns show redshifts of 2, 2.5, and 3. Physical coordinates 
are converted to comoving coordinates using the fiducial cosmology. 



The effective signal to noise ratio of individual observations spans about an order of magni- 
tude (see Figure 14) due to varying observing conditions, so it is important that observations 
are correctly weighted in the combination. However, an overall normalization error in the es- 
timated covariances C m still leads to a correctly weighted combination, even if the individual 
C m have different structures. 

3.2 Internal Covariance Tests and Refinement 

We exploit the large number of observations available to perform some internal cross checks 
on the individual C m and obtain improved estimates, as described below. First, we compare 
each observation d m to the combined observations using the chi-square statistic 

X 2 m = (dj - d)\C m - C)-\d j - d) (3.6) 
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Figure 13. Structure of the estimated covariance matrix C. The top-left panel (a) shows the ordered 
cigcnmodes X r of C in data (blue curve) and simulation (points). The top- right panel (b) shows 
the relative normalization of the 18 block-diagonal submatrices CV,-) as a function of separation AOj 
(black), as well as the relative normalization of the three 28 x 28 sub-matrices Cu^) corresponding to 
rcdshifts 2 (blue circles), 2.5 (green squares), and 3 (red diamonds), averaged over relative velocities 
Avi. The bottom-left panel (c) shows the common structure of the 18 C~ l blocks for each separation 
AOj. The bottom-right panel (d) shows the relative inverse covariance along the diagonal of this 
common block, as a function of relative velocity Avi/c, for redshifts 2 (blue circles), 2.5 (green 
squares), and 3 (red diamonds). 



where C m — C is the expected covariance of d m — d 

((d m - d)(d m - d)') = (d m d^) - 2(d m d<) + (dd<) 
= C m — 2C + ^ ^ CC m }C 

m' 

= C m -C, (3.7) 

and we have assumed that each d m is a Gaussian realization of some common true correlation 
function d(#o) with additional noise e m sampled from C m 

d m = d(6» ) + e m . (3.8) 

Figure 14 shows that, on average, our covariance estimates C m underestimate the actual 
variance seen in the data vectors d m by about 20%, but that this factor is independent of 
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Figure 14. Covariance studies of the 817 independent correlation function estimates in BOSS DR9. 
The top-left panel shows the distribution of covariance weights log 10 (| Cj\)/N for each observation 
j. The Xj values for each of the 817 observations are shown in the top-right panel as a function of 
observation index j, which is roughly ordered by date of observation, compared with their mean of 1.2. 
The bottom-left panel shows the distribution of normalized Xji defined in eqn. (3.6), for 1512 degrees 
of freedom, compared with their expected distribution (red curve) for Gaussian statistics with an 
overall normalization factor of 1.2 applied. The bottom-right panel shows the evolution of (x^) with 
eigenvalue rank r. Figure 15 repeats the bottom two panels after applying the rescaling procedure 
described in the text. 



when the observation was taken. We then search for a correlation between the degree of 
underestimation and the amount of predicted variance by diagonalizing each C m — C 

Cm — C = XmAmX^ (3-9) 

where A m is a diagonal matrix of the N = 1512 eigenvalues ranked in increasing size (larger 
predicted variance) 

{A-m)rs — &rs ^m.r j (3.10) 

and X m is a matrix whose columns are the corresponding eigenvectors of C m — C. We can 

.2 



then rewrite x? as 



where 

Xm,r = \ u m,r (3-12) 

and u m ^ r are the components of d m — d in the eigenvector basis for C m — C 

u m = X t m (d m - d) . (3.13) 
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Each Xm,r should be chi-square distributed with one degree of freedom (if our assumptions 
of Gaussianity are valid) so we use their mean over the M = 817 observations as a function 
of eigenvalue rank r 

r 

to provide another internal test, where we expect (Xr) — 1 for an eigenmode whose variance A r 
is correctly estimated. The results show (see Figure 14) that most of the 20% underestimation 
is due to the first ~ 250 eigenmodes with the largest predicted variances (but not including 
the 36 modes with artificially large variance mentioned above, which appear here with very 
small values of (Xr))- 

We can partially correct for the deviations with eigenvalue rank shown in Figure 14 
by independently rescaling each eigenmode r of C m by (x^) -1 - This procedure is not exact, 
however, since (Xr) _1 is derived from C m —C rather than C m (but we always have |C m | 3> |C|) 
and because each C m generally has different eigenvectors (but we expect similar eigenvectors 
since the C m differ primarily in their normalizations). The results of this rescaling procedure 
are shown in Figure 15 and demonstrate that all eigenmodes now have internally consistent 
d m fluctuations and covariance estimates C m . Two caveats to this procedure are that we do 
not rescale the 36 previously marginalized modes and that we impose the expected block- 
diagonal structure on the rescaled C m , eliminating any off-diagonal contributions introduced 
by numerical round-off errors. Note that the distribution of the Xm statistic has improved 
somewhat, compared with the Gaussian expectation, after rescaling. Finally, we combine the 
rescaled C m using eqn. (3.5) to obtain the final d and C used for the fits described below 
(and in Figs. 12 and 13). 




0.7 0.8 0.9 1.0 1.1 1.2 200 400 600 800 1000 1200 1400 



Observation ^,/d.o.f . Eigenmode Rank r 

Figure 15. Results of the rescaling procedure described in the text applied to BOSS DR9. Plots 
shows the same quantities as the bottom row of Figure 14 after applying the rescaling procedure 
described in the text. 

We use bootstrap sampling in two different ways in our analysis. First, we boostrap 
our 817 estimates of the data vector d m to obtain an estimate of their combined covariance 
that is independent of the C m , as a cross-check of our combination using eqn. (3.5). With 
M = 817 and 7Y = 1512, This approach is only feasible because of the 18 x 84 2 block-diagonal 
structure of C: we effectively bootstrap the 84(84 + l)/2 = 3570 elements of each block Ctj\ 
independently using 817 sub-vectors of length 84. We compare fit results obtained with 
the calculated combination C and the corresponding bootstrap estimate of C to check the 
validity of our calculated C. 
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The second type of bootstrap analysis we perform is to fit a large number of bootstrap 
samples and then compare the resulting distributions of best-fit parameter values with the 
parabolic errors calculated from the likelihood surface of the combined fit using d and C. To 
the extent that the errors in our fit parameters are not Gaussian, we expect some discrepancy 
and then prefer the bootstrap estimates. Since our observations span an order of magnitude 
in weights \C m \, a sampling with replacement generally yields a set of observations whose 
combined covariance does not match the combined covariance of all observations. Therefore, 
some care is required in assigning the total covariance used to fit each bootstrap sample since, 
although we ignore the resulting fit errors, we still require a correct relative weighting over 
the physical coordinate grid in eqn. (3.4). 

We generalize the combination method of eqn. (3.5) for bootstrap sampling to 

M M 
C 1 = ^ ] n m C m ^ , d = C 'y ^ n m C m ^ d m , (3.15) 

m=l m=l 

where n m > is the number of repetitions of observation m and X^m=i n m = W is the 
bootstrap sample size (usually M = M'). The resulting C is not the covariance of d since 
it incorrectly reduces the variance for double-counted observations. Instead, the correct 
covariance C to use is 

C = CD' 1 C (3.16) 

with 

M 

= n m C m ■ ( 3 - 17 ) 

m=l 

Note that eqn. (3.15) reduces to eqn. (3.5) when each = n m , so that each sample is either 
not used, n m = (if M' < M), or used exactly once, n m = 1. In the limit that all C m are 
identical, we find 

(X 2 ) M' + M-l 

W) ~ M ' (3 - 18) 

where x 2 is the chi-squared that would be obtained using C~ l instead of C~ l in eqn. (3.4). 
In the usual case of M' = M, we find a ratio 2 — 1/M ~ 2. 

Before performing fits, we generally apply some final cuts to select a subset of the iV 
three-dimensional grid points. Cuts can either be in physical or (using the fiducial cosmology) 
comoving coordinates. We implement cuts by first combining all observations using eqn. (3.5), 
without cuts, and then eliminating elements from the final data vector d, as well as the 
corresponding rows and columns of C. Note that the resulting C _1 is different from what 
we would have obtained by eliminating rows and columns from C _1 directly. Our nominal 
final cuts on physical coordinates are 0.003 < Av{/c < 0.083 and 5 < A6j < 165 arcminutes. 
We also cut on the co- moving separation 50 < r < 190 Mpc/h. After all final cuts, the size 
of our data vector is reduced from N = 3 x 504 = 1512 to 341 (z = 2) + 310 (z = 2.5) + 290 
{z = 3) = 941. 

4 Results 

In this section, we present results related to the fitting methods and what they reveal about 
the BOSS DR9 Lyman-a forest sample. Cosmological fit results derived from the same 
sample are presented in the companion paper ref. [29]. 
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4.1 Parameter Sensitivities 

In order to quantify and visualize how a dataset constrains each parameter p, it is instructive 
to plot the vector 

F p = d(0), p o^-M^) (4.1) 

where d(#) jP is the partial derivative of the theory prediction at some point in parameter space 
9, and o represents the Hadamard (entrywise) product. The sum of the resulting components 
is the Fisher information for parameter p [36] that specifies the expected maximum-likelihood 
error a p (in an ensemble-average sense) when all other parameter values are known 

^F Piljk ~a(p)- 2 . (4.2) 
ijk 

The theory derivatives d(0) tP encode the parameter sensitivity inherent to the model while 
the covariance C encodes the performance of a particular analysis method. The quantity F p 
is plotted for our baseline fit to BOSS DR9 for the isotropic scale parameter ai so in Figure 16 
and for the anisotropic scale parameters in Figures 17 (an) and 18 (a±). 




i iji iiii j ii^ i II 1 1 1 1 1 1 i ' i i i i i i 

U 50 100 150 °0 50 100 150 °0 50 100 150 

r ± (Mpc/h) rj_(Mpc/h) rj_(Mpc/h) 




50 100 150 50 100 150 50 100 150 



r±(Mpc/h) /-j^Mpc/h) r±(Mpc/h) 

Figure 16. Plots of the components of F p for the isotropic scale parameter a; so . The top row shows 
the intrinsic model sensitivity using C = |C| 1/n , while the bottom row uses the combined covariance 
matrix for BOSS DR9. The three plots on each row show element values in the (rii , r±) for z = 2 (left), 
2.5 (middle), and 3 (right), using the same color scale to preserve the relative magnitudes between 
redshifts. 

The top row of Figure 16 uses C = \C\ l l n to isolate the intrinsic model sensitivity, 
resulting in a Fisher error prediction of cr{a\ so ) = 0.017, with inverse-variance contributions 
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Figure 17. Plots of the components of F p for the anisotropic scale parameter ay. See caption to 
Figure 16 for details. 



of 12%, 28%, 60% from z = 2, 2.5, 3, respectively, due to the observed large increase in tracer 
bias b 2 (z) with redshift governed by the parameter 7^2. Information on the BAO scale comes 
predominantly from the forward direction due to the large observed value of the redshift- 
distortion parameter /3, and is concentrated on two rings in co- moving separation on either 
side of the BAO peak, where the theory is most sensitive to the peak position. When the 
combined covariance C is taken into account, the Fisher error prediction increases slightly 
to a(ct[ so ) = 0.018, and the region of peak sensitivity moves from fj, = 1 to fj, ~ 0.9, mostly 
due to the trend shown in Figure 13(b). Intermediate redshifts now dominate the statistical 
power due to the distribution of forest pixel redshifts, with inverse-variance contributions of 
29%, 61%, 10% from z = 2, 2.5, 3, respectively. 

When floating the scales parallel and perpendicular to the line of sight, the Fisher 
predicted errors are cr(a||) = 0.024 (0.021) and a(a±) = 0.050 (0.047) for the combined 
analysis covariance C (using C = |C| 1/n ). Note that our actual covariance actually reduces 
the expected error on cr(a±), relative to C = |C| 1//n , because of the decreasing errors with 
separation angle shown in Figure 13b. The peak predicted sensitivities are [i ~ 0.95 for a\\ 
and \i ~ 0.70 for a±. 

The results above were obtained with r" = r in eqn. (2.34), so that scale factors are 
applied only to the peak position and not the cosmological broadband. If, instead, we use 
r" = r' , the resulting Fisher components for an isotropic fit are shown in Figure 19, resulting 
in a smaller predicted error of a(ai so ) = 0.010 (0.013) for the combined analysis covariance 
C (using C = I CI 1 /™). Comparing with Figure 16, we conclude that a fit using coupled 
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Figure 18. Plots of the components of F p for the anisotropic scale parameter a±. See caption to 
Figure 16 for details. 



transforms obtains a smaller error (when all other parameter values are fixed) by using 
information outside of the peak region. However, this benefit is offset by a larger dependence 
on the broadband distortion model, and the final errors are comparable when marginalizing 
over distortion parameters. For fits to DR9, we prefer to localize the BAO scale measurement 
to the peak region and minimize our dependence on the distortion model. With larger data 
samples and improved analyses, we expect to have better control of distortion effects and 
this choice should be revisited. 

4.2 Data Reductions 

The 1512 fit inputs already represent a substantial reduction of data from the millions 
of Lyman-a forest pixel pairs from which they are derived, while preserving essentially all 
of the cosmologically relevant information. We can take this process one step further using 
the interpolated models of Section 2.4 to reduce down to estimates of multipoles of the 
correlation function or power spectrum. There are two main challenges in this process: first, 
the input £ijk has large correlated errors that will propagate through and, second, the input 
^ijk is expected to include some broadband distortion that prevents direct comparison of 
reduced data with theory. Therefore, we first test our data reductions on a mock dataset 
based on our fiducial cosmology and with the same covariance as data, but no distortion, 
and then compare these to the corresponding results obtained from data. 

We use three different types of model-independent fits for data reduction. The first fit 
estimates the r 2 - weighted correlation function multipoles at z = 2.4 as a linear interpolation 
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Figure 19. Plots of the components of F p for the isotropic scale parameter a; so , using r" = r' instead 
of r" = r in eqn. (2.34). See caption to Figure 16 for details. 



between 22 co-moving separations (5-140 in steps of 5, followed by 150, 160, and 180 Mpc/h) 
using eqn. (2.77), for a total reduction from 1512 to 66 values. The second fit estimates the 
k- weighted power spectrum multipoles at z = 2.4 as a set 20 of equally spaced band powers 
covering k = 0.03-0.33 h/Mpc with Ak = 0.015 h/Mpc, using eqn. (2.80) with 0-th order B- 
splines, for a total reduction from 1512 to 60 values. Our final model assumes that the P(k) 
multipoles are identical except for the normalization factor eqn. (2.78), which is expected to 
be a good approximation if broadband distortion is either small or else experiences similar 
redshift-space distortion as the linear theory, and results in 20 final values. In all of our data 
reduction fits, we assume (3{zq) = 1.4, 7^2 = 3.8, and 7/3 = 0, and our results are normalized 
using 6(2:0) = —0.183. We have made the results of these data reductions publicly available 
(see Appendix A). 

Figure 20 shows the ranked eigenvalues of the parameter covariance matrices for fits 
to data and the un-distorted mock dataset, and demonstrates that the data and mock have 
essentially the same parameter covariance structure. For the purposes of visualization, it 
is helpful to project out the largest eigenmodes of the parameter covariance matrix, which 
eliminates the largest sources of correlated errors but can also introduce a projection dis- 
tortion of the reduced data (in addition to any distortion from the analysis method). The 
optimum number of modes to project out is a balance between minimizing correlated errors 
and minimizing projection distortion effects. Figs. 21-23 show the results of projecting dif- 
ferent numbers of eigenvalues on the monopole of our un-distorted mock, and allows us to 
identify an appropriate number of modes to filter for each data reduction fit: -14, -36, and 
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Figure 20. Ranked eigenvalues of parameter covariance matrices for the three data reduction fits 
described in the text with, left to right, 66, 60, and 20 parameters, respectively. Red curves show- 
eigenvalues for the fits to data and points show fits to a mock dataset without broadband distortion. 
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Figure 21. Data reduction fit to an un-distorted mock with an increasing number of the largest 
cigcnmodes (shown in the top-left corner of each plot) projected out of the resulting parameter 
covariance matrix. Points show the r 2 -weighted monopole £o( r ? z) at z — 2.4 of the reduced data with 
errors calculated from the corresponding diagonal parameter covariance matrix elements. Curves 
show the fiducial model at z = 2.4 with (red) and without (blue, dashed) the projection distortion. 
Neither curve is a fit to the data points. 



Figs. 24-26 show the same reductions as Figs. 21-23, but applied to data instead of 
un-distorted mocks. Comparing these, we conclude that the projection distortion effects are 
similar for each reduction. We also note that the analysis introduces significant broadband 
distortion that is most apparent (especially with the r 2 weighting) as a broad excess above 
110 Mpc/h in the correlation-function monopole (Figure 21 with 14 modes projected out). 
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Figure 22. Data reduction fit to an un-distorted mock with an increasing number of the largest 
cigenmodes (shown in the top-left corner of each plot) projected out of the resulting parameter 
covariance matrix. Points show the fc-weighted deviation in the monopole AP(k, z) at z = 2.4 of the 
reduced data with errors calculated from the corresponding diagonal parameter covariance matrix 
elements. Curves show the fiducial model at z = 2.4 with (red) and without (blue, dashed) the 
projection distortion. Neither curve is a fit to the data points. 



Since this broadband distortion corresponds to \ow-k modes, the P{k) data reductions are 
relatively immune to it. Instead, they show an enhancement of the BAO oscillation signal 
relative to the fiducial model. Further study with additional data is needed to determine if 
this excess is a fortunate statistical fluctuation or if perhaps the signal is larger than expected. 

5 Discussion 

In this paper we have described a near-optimal method for fitting correlations observed in 
the Lyman-a forest of high-redshift quasars to measure the properties of the baryon acoustic 
oscillation feature. We also explore the properties of the BOSS DR9 correlation estimates 
described in ref. [29], and provide the code and input files necessary to reproduce our main 
results. 

Comparing with the Lyman-a fitting methods used in ref. [24], where correlations are 
modeled in the monopole and quadrupole at a single mean redshift, we have developed a 
fully three-dimensional fitting technique and included the effects of anisotropic non-linear 
broadening. These developments do not yield a significant improvement in the cosmological 
constraints that are possible with DR9, but provide more flexibility to study and quantify 
systematics and allow direct study the redshift evolution of parameters such as b(z) and @(z), 
and to accurately determine the effective redshift of our cosmological constraints. 
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Figure 23. Plots correspond to Figure 22 but use a data-reduction fit that assumes identical P(k) 
multipoles, resulting in 20 final parameters, as described in the text. 



The challenge of a three-dimensional analysis is that the resulting large covariance 
matrix is difficult to estimate and validate, especially with the limited number of mocks 
available for the DR9 Lyman-a forest. However, we have developed a novel method for 
internally testing and refining our large covariance matrix, as described in Section 3.2. The 
resulting covariance matrix still has large correlations that visually obscure the BAO signal 
but have little impact on our ability to measure its properties. We describe data reductions 
that retain most of the cosmologically relevant information and show how the expected 
features are visually revealed when large eigenmodes are projected out. 

In building models suitable for a three-dimensional analysis, we find that care is needed 
to isolate the BAO feature in the correlation function and describe an alternative to the 
"no- wiggles" approach of ref. [34] that is better suited for this purpose. When considering 
the line of sight (aii) and transverse (a±) scale factors independently, we find that a first- 
order treatment is numerically inaccurate near the BAO peak, and identify the important 
second-order terms. We introduce anisotropic non-linear broadening using an approximation 
that is valid for a large range of (3. We also use a highly flexible parametrization of both 
multiplicative and additive broadband distortions to accomodate the expected systematics 
of the correlation- function estimates due to effects such as a biased continuum estimate. 
Finally, we demonstrate that information on the BAO feature in the Lyman-a forest arises 
primarily along the line of sight, as expected due to the large redshift-space distortion, but 
that the contributions from H(z) and Da{z) can still be independently measured. 

The first measurements of baryon acoustic oscillations in the Lyman-a forest [24, 29] 
provide a new confirmation of our basic picture of cosmological evolution dominated by dark 
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Figure 24. Plots correspond to Figure 21 but are based on data instead of the un-distorted mock. 



energy, but do not yet significantly constrain the parameters of the simplest cosmological 
models. With the SDSS-III high-redshift quasar sample expected to triple the DR9 statistics 
over the next 18 months, that situation is likely to change and will require near-optimal 
fitting methods and a flexible framework for careful study of systematics. 
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Figure 25. Plots correspond to Figure 22 but are based on data instead of the un-distorted mock. 
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Figure 26. Plots correspond to Figure 23 but are based on data instead of the un-distorted mock. 
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A Public Access to Data and Code 

The software used to generate the results in this paper and ref. [29] are publicly available at 
http://github.com/baofit/. We also provide instructions to install and run the software, 
together with the BOSS DR9 correlation-function estimates and configuration files necessary 
to reproduce our main results, at http://darkmatter.ps.uci.edu/baofit/. The software 
is written in C++ and uses MINUIT [37] for likelihood minimization. Data and configuration 
files are in plain text format. 
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